function [vk_dH, vep_dH, ck_dH, cu_dH] ...
    = fn_solve_tba_dH(T, H, vk, vep, ck, cu)


Num_point = 400;

Nc = length(ck);
cP = cell(1, Nc);
cQ = cell(1, Nc);

for n = 1:Nc
    
    Pn = fn_u2Pn(n, H, T, cu{n});
    cP{n} = Pn;
    Qn = fn_u2Qn(n, H, T, cu{n});
    cQ{n} = Qn;
    
end

ck_dH = cell(1, Nc);
cu_dH = cell(1, Nc);

for n = 1:Nc
    
    ck_dH{n} = linspace(0, 1000, 200);
    cu_dH{n} = zeros(1, length(ck_dH{n}));
    
end

vk_dH = vk;
vep_dH = -1/2*tanh(H/2/T)*ones(1, length(vk));

while 1
    
    vep_dH0 = vep_dH;
    ck_dH0 = ck_dH;
    cu_dH0 = cu_dH;
    
    tvep_dH = zeros(1, length(vk));
    for i = 1:length(vk_dH)
        
        k = vk_dH(i);
        res = -1/2*tanh(H/2/T);
%         res = -1;
        res = res + fn_int_an(2, vk, fn_G(vep, T).*vep_dH, k);
            
        for n = 1:Nc
            
            [tx2, tf2] = fn_combine(ck{n}, cP{n},  cQ{n}, ck_dH{n}, cu_dH{n});
            res = res - T*fn_int_an(n, tx2, tf2, k);
            
        end
        
        tvep_dH(i) = res;
        
    end
    
    vep_dH = tvep_dH;
    
    vF = -1/T*fn_G(vep, T).*vep_dH;
    
    vnorm = zeros(1, Nc);
    for n = 1:Nc
        
        [vp, vf] ...
            = fn_iteration2(n, vk, vF, ck, cP, cQ, ck_dH, cu_dH, Num_point);
        ck_dH{n} = vp;
        cu_dH{n} = vf;
        
        vnorm(n) = fn_norm_2(ck_dH{n},cu_dH{n},ck_dH0{n},cu_dH0{n});
        
    end
    
    diss1 = fn_norm(vep_dH - vep_dH0);
    diss2 = fn_norm(vnorm);
    
    fprintf('diss = %1.2e   ', diss1);
    fprintf('diss = %1.2e   \n', diss2);
    
    if diss1<1e-8 && diss2<1e-8
        break;
    end
    
end

end